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Abstract 

Many processes in biology, from the regulation of gene expression in 
bacteria to memory in the brain, involve switches constructed from net- 
works of biochemical reactions. Crucial molecules are present in small 
numbers, raising questions about noise and stability. Analysis of noise 
in simple reaction schemes indicates that switches stable for years and 
switchable in milliseconds can be built from fewer than one hundred 
molecules. Prospects for direct tests of this prediction, as well as im- 
plications, are discussed. 



1 Introduction 

The problem of building a reliable switch arises in several different biological 
contexts. The classical example is the switching on and off of gene expression 
during development jj], , or in simpler systems such as phage A || || . It is likely 
that the cell cycle should also be viewed as a sequence of switching events among 
discrete states, rather than as a continuously running clock ||. The stable 
switching of a specific class of kinase molecules between active and inactive 
states is believed to play a role in synaptic plasticity, and by implication in the 
maintenance of stored memories |g, 0]. Although many details of mechanism 
remain to be discovered, these systems seem to have several common features. 
First, the stable states of the switches are dissipative, so that they reflect a 
balance among competing biochemical reactions. Second, the total number of 
molecules involved in the construction of the switch is not large. Finally, the 
switch, once flipped, must be stable for a time long compared to the switching 
time, perhaps — for development and for memory — even for a time comparable 
to the life of the organism. Intuitively we might expect that systems with small 
numbers of molecules would be subject to noise and instability S, and while 
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this is true we shall see that extremely stable biochemical switches can in fact 
be built from a few tens of molecules. This has interesting implications for how 
we think about several cellular processes, and should be testable directly. 

Many biological molecules can exist in multiple states, and biochemical 
switches use this molecular multistability so that the state of the switch can 
be 'read out' by sampling the states (or enzymatic activities) of individual 
molecules. Nonetheless, these biochemical switches are based on a network of 
reactions, with stable states that are collective properties of the network dynam- 
ics and not of any individual molecule. Most previous work on the properties 
of biochemical reaction networks has involved detailed simulation of particular 
kinetic schemes ||, [lO), for example in discussing the kinase switch that is in- 
volved in synaptic plasticity |pd[| . Even the problem of noise has been discussed 



hcuristically in this context 1 12 . The goal in the present analysis is to separate 
the problem of noise and stability from other issues, and to see if it is possible 
to make some general statements about the limits to stability in switches built 
from a small number of molecules. This effort should be seen as being in the 
same spirit as recent work on bacterial chemotaxis, where the goal was to un- 
derstand how certain features of the computations involved in signal processing 
can emerge robustly from the network of biochemical reactions, independent of 
kinetic details fbl. 



2 Stochastic kinetic equations 

Imagine that we write down the kinetic equations for some set of biochemical 
reactions which describe the putative switch. Now let us assume that most of 
the reactions are fast, so that there is a single molecular species whose con- 
centration varies more slowly than all the others. Then the dynamics of the 
switch essentially are one dimensional, and this simplification allows a complete 
discussion using standard analytical methods. In particular, in this limit there 
are general bounds on the stability of switches, and these bounds are indepen- 
dent of (incompletely known) details in the biochemical kinetics. It should be 
possible to make progress on multidimensional versions of the problem, but the 
point here is to show that there exists a limit in which stable switches can be 
built from small numbers of molecules. 

Let the number of molecules of the 'slow species' be n. All the different 
reactions can be broken into two classes: the synthesis of the slow species at 
a rate f(n) molecules per second, and its degradation at a rate g{n) molecules 
per second; the dependencies on n can be complicated because they include the 
effects of all other species in the system. Then, if we could neglect fluctuations, 
we would write the effective kinetic equation 

f£ = /(n)- fl (n). (1) 
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If the system is to function as a switch, then the stationarity condition f(n) = 
g(n) must have multiple solutions with appropriate local stability properties. 

The fact that molecules are discrete units means that we need to give the 
chemical kinetic Eq. (1) another interpretation. It is the mean field approxima- 
tion to a stochastic process in which there is a probability per unit time f(n) 
of making the transition n — > n + 1, and a probability per unit time g(n) of 
the opposite transition n — ► n — 1. Thus if we consider the probability P(n,t) 
for there being n molecules at time t, this distribution obeys the evolution (or 
'master') equation 

^^ = /(n-l)P(n-l,t)+ S (n+l)P(n+l,t)-[/(n)+ fl (n)]P(n,*) ) (2) 

with obvious corrections for n = 0, 1. We are interested in the effects of stochas- 
ticity for n not too small. Then 1 is small compared with typical values of n, and 
we can approximate P(n,t) as being a smooth function of n. We can expand 
Eq. (Q) in derivatives of the distribution, and keep the leading terms: 

= -L { [9{n) " f ^ P ^ + l^Un+9(n)}P(n,t^ . (3) 

This is analogous to the diffusion equation for a particle moving in a potential, 
but this analogy works only if allow the effective temperature to vary with the 
position of the particle. 

As with diffusion or Brownian motion, there is an alternative to the diffusion 
equation for P(n, t) and this is to write an equation of motion for n(t) which 
supplements Eq. (jlj) by the addition of a random or Langevin force 

J = f(n)-g(n)+m, (4) 
mm) = [f(n) + g(n)]S(t-t'). (5) 

From the Langevin equation we can also develop the distribution functional 
for the probability of trajectories n(t). It should be emphasized that all of 
these approaches are equivalent provided that we are careful to treat the spatial 
variations of the effective temperature |Q .[] In one dimension this complication 
does not impede solving the problem. For any particular kinetic scheme we 
can compute the effective potential and temperature, and kinetic schemes with 
multiple stable states correspond to potential functions with multiple minima. 

1 In a review written for a biological audience, McAdams and Arkin |Ie| ] state that Langevin 
methods are unsound and can yield invalid predictions precisely for the case of bistable reaction 
systems which interests us here; this is part of their argument for the necessity of stochastic 
simulation method s a s opposed to analytic approaches. Their reference for the failure of 
Langevin methods 1161 , however, seems to consider only Langevin terms with constant spectral 
density, thus ignoring (in the present language) the spatial variations of effective temperature. 
For the present problem this would mean replacing the noise correlation function [/(n) + 
g(n)]S(t — t') in Eq. (H) by Q5(t — t') where Q is a constant. This indeed is wrong, and is not 
equivalent to the master equation. On the other hand, if the arguments of Refs. |l6| were 
generally correct, they would imply that Langevin methods could not used for the description 
of Brownian motion with a spatially varying temperature, and this would be quite a surprise. 
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3 Noise induced switching rates 



We want to know how the noise term destabilizes the distinct stable states of 
the switch. If the noise is small, then by analogy with thermal noise we ex- 
pect that there will be some small jitter around the stable states, but also some 
rate of spontaneous jumping between the states, analogous to thermal activa- 
tion over an energy barrier as in a chemical reaction. This jumping rate should 
be the product of an "attempt frequency" — of order the relaxation rate in the 
neighborhood of one stable state — and a "Boltzmann factor" that expresses the 
exponentially small probability of going over the barrier. For ordinary chemical 
reactions this Boltzmann factor is just exp(— F' /ksT), where FT is the activa- 
tion free energy. If we want to build a switch that can be stable for a time much 
longer than the switching time itself, then the Boltzmann factor has to provide 
this large ratio of time scales. 

There are several ways to calculate the analog of the Boltzmann factor for 
the dynamics in Eq. (Q). The first step is to make more explicit the analogy 
with Brownian motion and thermal activation. Recall that Brownian motion of 
an overdamped particle is described by the Langevin equation 

7§=-V"dO+T?(t), (6) 

where 7 is drag coefficient of the particle, V(x) is the potential, and the noise 
force has correlations (r)(t)r](t')) = 2-fTS(t — t'), where T is the absolute tem- 
perature measured in energy units so that Boltzmann's constant is equal to one. 
Comparing with Eq. (0), we see that our problem is equivalent to a particle 
with 7 = 1 in an effective potential V e g(n) such that V^ s {n) = g(n) — f(n), at 
an effective temperature T c tf(n) = [f(n) + g(n)]/2. 

If the temperature were uniform then the equilibrium distribution of n would 
be P eq (n) oc exp[— V c e {n)/T c s\. With nonuniform temperature the result is (up 
to weakly varying prcfactors) 

P cq (n) oc exp[-U(ri)] (7) 

"<»> = £ d * v M _ (8) 

One way to identify the Boltzmann factor for spontaneous switching is then to 
compute the relative equilibrium occupancy of the stable states (n and nx) and 
the unstable "transition state" at n* . The result is that the effective activation 
energy for transitions from a stable state at n = tiq to the stable state at 
n = ri\ > no is 

FHno - n x ) = 2k B T P dn g [ w) -ff"j , (9) 
J no 9{n) + /(«) 

where n* is the unstable point, and similarly for the reverse transition, 

Ftfa - n ) = 2k B T I" d /W-ffN (10) 
Jn t 9{n) + /(n) 
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An alternative approach is to note that the distribution of trajectories n(t) 
includes locally optimal paths that carry the system from each stable point up 
to the transition state; the effective activation free energy can then be written as 
an integral along these optimal paths. The use of optimal path ideas in chemical 
kinetics has a long history, going back at least to Onsager. A discussion in the 
spirit of the present one is Ref . O] . For equations of the general form 

^ = -y e ' ff (n)+£(i), (11) 

with (£(£)£(£')) = %Tes(t)5(t — t'), the probability distribution for trajectories 
P[n(t)] can be written as (m| 

P[n{t)] = exp(-S[n(f)]) (12) 

S[n(t)} = \f dtjr^{h(t) + V: s (n(t))} 2 -If dtV^(n(t)). (13) 

If the temperature T e ff is small, then the trajectories that minimize the action 
should be determined primarily by minimizing the first term in Eq. ([l3|), which 
is ~ 1/T c ff. Identifying the effective potential and temperature as above, the 
relevant term is 

1 f dt [n-f(n) +9 (n) ? = W n> 1 f JJjn) - 9 (n)f 



2 J f(n)+g(n) 2 J f(n)+g(n) 2 J f(n)+g(n) 

dtn — — . 14 

f{n) + g[n) 

We are searching for trajectories which take n(t) from a stable point no where 
/(no) = g(no) through the unstable point where / and g are again equal but 
the derivative of their difference (the curvature of the potential) has changed 
sign. For a discussion of the analogous quantum mechanical problem of tunnel- 
ing in a double well, see Ref. [^8). First we note that along any trajectory from 
no to n* we can simplify the third term in Eq. (|l4|): 

, . f(n) — q(n) f n * , f(n) — q(n) , . 

dth J ) ' y )[ = / dn J ) : ' y ) . 15 
f{n)+g(n) J no f(n)+g{n) 

This term thus depends on the endpoints of the trajectory and not on the path, 
and therefore cannot contribute to the structure of the optimal path. In the 
analogy to mechanics, the first two terms are equivalent to the (Euclidean) 
action for a particle with position dependent mass in a potential; this means 
that along extremal trajectories there is a conserved energy 

E= l l [/(n)-. 9 (n)] 2 

2f(n)+g(n) 2 f(n)+g(n) ' ( ' 

At the endpoints of the trajectory we have n = and /(n) = g(n), and so we 
are looking for zero energy trajectories, along which 

n(t) = ±[f(n(t)) - g(n(t))}- (17) 
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Substituting back into Eq. (pjj), and being careful about the signs, we find once 
again Eq's. (|,|l^). 

Both the 'transition state' and the optimal path method involve approxi- 
mations, but if the noise is not too large the approximations are good and the 
results of the two methods agree. Yet another approach is to solve the master 
equation (||) directly, and again one gets the same answer for the switching rate 
when the noise is small, as expected since all the different approaches are all 
equivalent if we make consistent approximations. It is much more work to find 
the prefactors of the rates, but we are concerned here with orders of magnitude, 
and hence the prefactors aren't so important. 



4 Interpretation 

The crucial thing to notice in this calculation is that the integrands in Eq's. 
(p|JlC|) are bounded by one, so the activation energy (in units of the thermal 
energy ksT) is bounded by twice the change in the number of molecules. Trans- 
lating back to the spontaneous switching rates, the result is that the noise driven 
switching time is longer than the relaxation time after switching by a factor that 
is bounded, 

spontaneous switching time , A . ,„ . 

-i = — -2 < exp An , 18 

relaxation time 

where An is the change in the number of molecules required to go from one 
stable 'switched' state to the other. Imagine that we have a reaction scheme 
in which the difference between the two stable states corresponds to roughly 25 
molecules. Then it is possible to have a Boltzmann factor of up to exp(25) ~ 
10 10 . Usually we think of this as a limit to stability: with 25 molecules we 
can have a Boltzmann factor of no more than ~ 10 10 . But here I want to 
emphasize the positive statement that there exist kinetic schemes in which just 
25 molecules would be sufficient to have this level of stability. This corresponds 
to years per millisecond: with twenty five molecules, a biochemical switch that 
can flip in milliseconds can be stable for years. Real chemical reaction schemes 
will not saturate this bound, but certainly such stability is possible with roughly 
100 molecules. The genetic switch in A phage operates with roughly 100 copies 
of the repressor molecules (^], and even in this simple system there is extreme 
stability: the genetic switch is flipped spontaneously only once in 10 5 generations 
of the host bacterium |Q . Kinetic schemes with greater cooperativity get closer 
to the bound, achieving greater stability for the same number of molecules. 

In electronics, the construction of digital elements provides insulation against 
fluctuations on a microscopic scale and allows a separation between the logical 
and physical design of a large system. We see that, once a cell has access to sev- 
eral tens of molecules, it is possible to construct 'digital' switch elements with 
dynamics that are no longer significantly affected by microscopic fluctuations. 
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Furthermore, weak interactions of these molecules with other cellular compo- 
nents cannot change the basic 'states' of the switch, although these interactions 
can couple state changes to other events. 

The importance of this 'digitization' on the scale of 10 — 100 molecules is 
illustrated by different models for pattern formation in development. In the 
classical model due to Turing, patterns are expressed by spatial variations in 
the concentration of different molecules, and patterns arise because uniform 
concentrations are rendered unstable through the combination of nonlinearities 
in the kinetics with the different diffusion constants of different substances. In 
this picture, the spatial structure of the pattern is linked directly to physical 
properties of the molecules. An alternative that each spatial location is labelled 
by a set of discrete possible states, and patterns evolve out of the 'automaton' 
rules by which each location changes state in relation to the neighboring states. 
In this picture states and rules are more abstract, and the dynamics of pattern 
formation is really at a different level of description from the molecular dynamics 
of chemical reactions and diffusion. Reliable implementations of automaton 
rules apparently are accessible as soon as the relevant chemical reactions involve 
a few dozen molecules. 

Biochemical switches have been reconstituted in vitro, but I am not aware 
of any attempts to verify that stable switching is possible with small numbers 
of molecules. It would be most interesting to study model systems in which one 
could confine and monitor sufficiently few molecules that it becomes possible 
to observe spontaneous switching, that is the breakdown of stability. Although 
genetic switches have certain advantages, even the simplest systems would re- 
quire full enzymatic apparatus for gene expression (but see Ref. plj for recent 
progress on controllable in vitro expression systems).^] Kinase switches are much 
simpler, since they can be constructed from just a few proteins and can be trig- 
gered by calcium; caged calcium allows for an optical pulse to serve as input. 

At reasonable protein concentrations, 10 — 100 molecules are found in a 
volume of roughly 1 (fim) 3 . Thus it should be possible to fabricate an array of 
'cells' with linear dimensions ranging from 100 nm to 10 /jm, such that solutions 
of kinase and accessory proteins would switch stably in the larger cells but 
exhibit instability and spontaneous switching in the smaller cells. The state 

2 Note also that reactions involving polymer synthesis (mRNA from DNA or protein from 
mRNA) are not 'elementary' reactions in the sense described by Eq. (fe|). Synthesis of a 
single mRNA molecule involves thousands of steps, each of which occurs (conditionally) at 
constant probability per unit time, and so the noise in the overall synthesis reaction is very 
different. If the synthesis enzymes are highly processive, so that the polymerization apparatus 
incoporates many monomers into the polymer before 'backing up' or falling off the template, 
then synthesis itself involves a delay but relatively little noise; the dominant source of noise 
becomes the assembly and disassembly of the polymerization complex. Thus there is some 
subtlety in trying to relate a simple model to the complex sequence of reactions involved 
in gene expression. On the other hand a detailed simulation is problematic, since there are 
so many different elementary steps with unknown rates. This combination of circumstances 
would make experiments on a minimal, in vitro genetic switch espcially interesting. 
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of the switch could be read out by including marker proteins that would serve 
as substrates of the kinase but have, for example, fluorescence lines that are 
shifted by phosphorylation, or by having fluorescent probes on the kinase itself; 
transitions of single enzyme molecules should be observable |l|, |2Cj . 

A related idea would be to construct vesicles containing ligand gate ion 
channels which can conduct calcium, and then have inside the vesicle enzymes 
for synthesis and degradation of the ligand which are calcium sensitive. The 
cGMP channels of rod photoreceptors are an example, and in rods the cy- 
clase synthesizing cGMP is calcium sensitive, but the sign is wrong to make a 



switch 1 22 ; presumably this could solved by appropriate mixing and matching 
of protein components from different cells. In such a vesicle the different stable 
states would be distinguished by different levels of internal calcium (as with 
adaptation states in the rod), and these could be read out optically using cal- 
cium indicators; caged calcium would again provide an optical input to flip the 
switch. Amusingly, a close packed array of such vesicles with ^100 nm dimen- 
sion would provide an optically addressable and writable memory with storage 
density comparable to current RAM, albeit with much slower switching. 

In summary, it should be possible to build stable biochemical switches from 
a few tens of molecules, and it seems likely that nature makes use of these. To 
test our understanding of stability we have to construct systems which cross the 
threshold for observable instabilities, and this seems accessible experimentally 
in several systems. 
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